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Convergence  Results  for  Ant  Routing  Algorithms 
via  Stochastic  Approximation  and  Optimization 


Punyaslok  Purkayastha  and  John  S.  Baras 


Abstract — “Ant  algorithms”  have  been  proposed  to  solve  a 
variety  of  problems  arising  in  optimization  and  distributed 
control.  They  form  a  subset  of  the  larger  class  of  “Swarm 
Intelligence”  algorithms.  The  central  idea  is  that  a  ’’swarm” 
of  relatively  simple  agents  can  interact  through  simple  mecha¬ 
nisms  and  collectively  solve  complex  problems.  Instances  that 
exemplify  the  above  idea  abound  in  nature.  The  abilities  of  ant 
colonies  to  collectively  accomplish  complex  tasks  have  served 
as  sources  of  inspiration  for  the  design  of  “Ant  algorithms”. 
Examples  of  “Ant  algorithms”  are  the  set  of  “Ant  Routing”  al¬ 
gorithms  that  have  been  proposed  for  communication  networks. 
We  analyze  in  this  paper  Ant  Routing  Algorithms  for  packet- 
switched  wireline  networks.  The  algorithm  retains  most  of  the 
salient  and  attractive  features  of  Ant  Routing  Algorithms.  The 
scheme  is  a  multiple  path  probabilistic  routing  scheme,  that  is 
fully  adaptive  and  distributed.  Using  methods  from  adaptive 
algorithms  and  stochastic  approximation,  we  show  that  the 
evolution  of  the  link  delay  estimates  can  be  closely  tracked  by 
a  deterministic  ODE  system.  A  study  of  the  equilibrium  points 
of  the  ODE  gives  us  the  equilibrium  behavior  of  the  routing 
algorithm,  in  particular,  the  equilibrium  routing  probabilities, 
and  mean  delays  in  the  links  under  equilibrium.  We  also  show 
that  the  fixed-point  equations  that  the  equilibrium  routing 
probabilities  satisfy  are  actually  the  necessary  and  sufficient 
conditions  of  an  appropriate  optimization  problem.  Simulations 
supporting  the  analytical  results  are  provided. 

I.  INTRODUCTION 

“Ant  algorithms”  constitute  a  class  of  algorithms  that  have 
been  proposed  to  solve  a  variety  of  problems  arising  in 
optimization  and  distributed  control.  They  form  a  subset 
of  the  larger  class  of  what  are  referred  to  as  “Swarm 
Intelligence”  algorithms,  a  topic  extensively  discussed  in  the 
book  by  Bonabeau,  Theraulez,  and  Dorigo  [5],  The  central 
idea  here  is  that  a  “swarm”  of  relatively  simple  agents  can 
interact  through  simple  mechanisms  and  collectively  solve 
complex  problems.  Instances  that  exemplify  the  above  idea 
abound  in  nature.  Bonabeau,  Theraulez,  and  Dorigo  [5] 
give  examples  of  insect  societies  like  those  of  ants,  honey 
bees,  and  wasps,  which  accomplish  fairly  complex  tasks  of 
building  intricate  nests,  finding  food,  responding  to  external 
threats  etc.,  even  though  the  individual  insects  themselves 
have  limited  capabilities.  The  abilities  of  ant  colonies  to 
collectively  accomplish  complex  tasks  have  served  as  sources 
of  inspiration  for  the  design  of  “Ant  algorithms”. 

Examples  of  “Ant  algorithms”  are  the  set  of  “Ant-Based 
Routing”  algorithms  that  have  been  proposed  for  commu- 
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nication  networks.  It  has  been  observed  in  an  experiment 
conducted  by  biologists  called  the  double  bridge  experiment 
[7],  that  a  colony  of  ants,  when  presented  with  two  paths 
to  a  source  of  food,  is  able  to  collectively  converge  to  the 
shorter  path.  Every  ant  lays  a  trail  of  a  chemical  substance 
called  pheromone  as  it  walks  along  a  path;  subsequent  ants 
follow  and  reinforce  this  trail.  This  leads  progressively  to  a 
large  accumulation  of  pheromone  on  the  shorter  path,  which 
is  how  ants  discover  the  shorter  path.  Most  of  the  Ant- 
Based  Routing  Algorithms  (called  Ant  Routing  Algorithms, 
for  short)  proposed  in  the  literature  are  inspired  by  this 
basic  idea.  These  algorithms  employ  probe  packets  called 
ant  packets  (analogues  of  ants)  to  explore  the  network 
and  measure  various  quantities  related  to  network  routing 
performance  like  link  and  path  delays.  These  measurements 
are  used  to  construct  and  update  the  routing  tables  at  the 
network  nodes.  The  update  algorithms  tend  to  reinforce  those 
outgoing  links  which  lead  to  paths  with  lower  delays. 

Schoonderwoerd  et.  al.  [11],  [5]  tested  an  Ant  Routing 
Algorithm  on  the  British  Telecom  telephone  network,  and 
reported  superior  performance  compared  to  other  algorithms 
including  shortest  paths  based  schemes.  This  generated  in¬ 
terest  in  the  study  of  Ant  Routing  Algorithms  for  both 
connection-oriented  networks  and  for  packet-switched  net¬ 
works  (both  wired  and  wireless).  Ant  Routing  Algorithms 
for  wireless  networks  have  been  proposed  and  analyzed  in 
papers  by  Baras  and  Mehta  [1]  and  Gunes  et.  al.  [9], 

We  consider  in  this  paper  packet-switched  wired  networks. 
Ant  Routing  Algorithms  for  such  networks  have  been  pro¬ 
posed  in  the  works  of  Di  Caro  and  Dorigo  [8]  and  Bean  and 
Costa  [2].  Though  a  large  number  of  ant  routing  algorithms 
have  been  proposed,  very  few  analytical  studies  are  available 
in  the  literature.  Examples  of  such  analytical  studies  can  be 
found  in  [12],  [10],  and  [2].  Yoo,  La,  and  Makowski  [12]  and 
Gutjahr  [10]  consider  very  simple  network  cases  where  the 
delays  in  the  links  of  the  network  are  deterministic  quantities, 
independent  of  the  offered  traffic  loads.  They  show  that 
ant  algorithms  converge  to  the  shortest  path  solutions  in 
the  network.  Bean  and  Costa  [2]  propose  a  multiple-path 
routing  scheme  which  tries  to  form  estimates  of  the  delays 
and  use  them  to  update  routing  tables.  The  link  delays 
can  be  stochastically  varying.  The  authors  employ  a  time- 
scale  separation  approximation  whereby  the  delay  estimates 
are  computed  before  the  routing  probabilities  are  updated. 
They  find  that  results  from  a  numerical  study  based  on  their 
model  and  those  from  simulations  agree  well.  However,  the 
exact  nature  of  the  time-scale  separation  is  not  clear  nor 
is  any  convergence  result  provided.  Borkar  and  Kumar  [6] 


use  stochastic  approximation  methods  to  prove  convergence 
of  their  scheme  called  Wardrop  routing.  Their  framework  is 
similar  to  ours  -  they  have  a  delay  estimation  scheme  and  a 
probability  update  scheme  which  utilizes  the  delay  estimates. 
Their  probability  update  scheme  moves  on  a  slower  time- 
scale  than  their  delay  estimation  scheme.  By  two  time-scale 
stochastic  approximation  methods  they  study  convergence  of 
their  algorithm. 

We  consider  the  algorithm  proposed  by  Bean  and  Costa 
[2]  in  this  paper.  The  algorithm  retains  most  of  the  attractive 
features  of  Ant  Routing  Algorithms.  The  scheme  is  fully 
adaptive  and  distributed.  We  consider  in  this  paper  the  simple 
routing  scenario  where  data  traffic  entering  a  single  source 
node  has  to  be  routed  to  a  single  destination  node,  and  there 
are  N  available  parallel  paths  between  them.  We  model  the 
arrival  processes  and  packet  lengths  of  both  the  ant  and  the 
data  streams  that  arrive  at  the  source  node,  and  argue,  using 
methods  from  adaptive  algorithms  and  stochastic  approxi¬ 
mation,  that  the  evolution  of  the  link  delay  estimates  can 
be  closely  tracked  by  a  deterministic  ODE  system,  when 
the  step  size  of  the  estimation  scheme  is  small.  Then  a 
study  of  the  equilibrium  points  of  the  ODE  gives  us  the 
equilibrium  behavior  of  the  routing  algorithm;  in  particular, 
the  equilibrium  routing  probabilities,  and  mean  delays  in  the 
N  paths  under  equilibrium  can  be  obtained.  We  also  show 
that  the  fixed-point  equations  that  the  equilibrium  routing 
probabilities  satisfy  are  actually  the  necessary  and  sufficient 
conditions  of  a  convex  minimization  problem. 

Our  paper  is  organized  as  follows.  Section  II  provides  the 
general  framework  of  ant  routing  and  the  routing  scheme  we 
consider.  Section  III  provides  a  detailed  description  of  our  N 
parallel  paths  model.  Section  IV  contains  an  analysis  of  the 
routing  scheme.  Section  V  provides  illustrative  simulation 
results,  and  in  Section  VI  we  provide  a  few  conclusive 
remarks. 

II.  GENERAL  FRAMEWORK  OF  ANT  ROUTING 
ALGORITHMS  AND  OUR  SCHEME 

We  provide,  in  this  section,  a  brief  description  of  the  gen¬ 
eral  framework  of  ant  routing  for  a  (wired)  communication 
network.  The  framework  that  we  follow  is  the  one  described 
in  Di  Caro  and  Dorigo  [7],  [8].  Alongside,  we  describe  the 
scheme  due  to  Bean  and  Costa  [2],  that  we  analyse  in  this 
paper. 

Every  node  i  in  the  network  maintains  two  key  data 
structures  -  a  matrix  of  routing  probabilities,  the  routing 
table  and  a  matrix  of  various  kinds  of  statistics,  called 
the  local  network  information  table,  C(i).  For  a  particular 
node  i,  let  k )  denote  the  set  of  neighbors  of  i  through 
which  node  i  routes  packets  towards  destination  />::.  For  a 
communication  network  consisting  of  M  nodes,  the  matrix 
of  routing  probabilities,  lZ(i ),  has  M  —  1  columns,  corre¬ 
sponding  to  the  M  —  1  destinations  towards  which  node  i 
could  route  packets,  and  M  —  1  rows,  corresponding  to  the 
maximum  number  of  neighbor  nodes  through  which  node  i 
could  route  packets  to  a  particular  destination.  The  entries 
of  lZ(i)  are  the  probabilities  denotes  the  probability 


of  routing  an  incoming  packet  at  node  i  and  bound  for 
destination  k  via  the  neighbor  j  €  A/"(z,  k).  The  matrix 
C{i)  has  the  same  dimensions  as  and  its  (),  j)-th  entry 
contains  various  statistics  pertaining  to  the  route  (i,  j, . . . ,  k). 
Examples  of  such  statistics  could  be  mean  delay  and  delay 
variance  estimates  of  the  route  ( i,j,...,k ).  These  statistics 
are  maintained  and  updated  based  on  the  information  the  ant 
packets  collect  about  the  route.  The  matrix  C(i)  represents 
the  characteristics  of  the  network  that  are  learned  by  the 
nodes  through  the  ant  packets,  based  on  which  local  decision¬ 
making,  and  the  updating  of  the  routing  table  IZ(i),  are  done. 
The  iterative  algorithms  that  are  used  to  update  C(i)  and 
IZ(i)  will  be  referred  to  as  the  learning  algorithms. 

We  now  describe  the  mechanism  of  operation  of  ant- 
based  routing  algorithms.  For  ease  of  exposition,  we  restrict 
attention  to  a  particular  fixed  destination  node,  and  consider 
the  problem  of  routing  from  every  other  node  to  this  node, 
which  we  label  as  D. 

Forward  ant  generation  and  routing.  At  certain  inter¬ 
vals,  forward  ant  (FA)  packets  are  launched  from  a  node  i 
towards  the  destination  node  D  to  discover  low  delay  paths  to 
it.  The  FA  packets  sample  walks  on  the  graph  representing 
the  communication  network,  based  on  the  current  routing 
probabilities  at  the  nodes.  FA  packets  share  the  same  queues 
as  data  packets  and  so  experience  similar  delay  characteris¬ 
tics  as  data  packets.  Every  FA  packet  maintains  a  stack  of 
data  structures  containing  the  IDs  of  nodes  in  its  path  and  the 
per  hop  delays  (or  other  relevant  information)  encountered. 
The  per  hop  delay  measurements  are  obtained  through  time 
stamping  of  the  packets  as  they  pass  from  the  various  nodes. 

Backward  ant  generation  and  routing.  Upon  arrival  of 
an  FA  at  the  destination  node  D ,  a  backward  ant  (BA)  packet 
is  generated.  The  FA  packet  transfers  its  stack  to  the  BA.  The 
BA  packet  then  retraces  back  to  the  source  the  path  traversed 
by  the  FA  packet.  BA  packets  travel  back  in  high  priority 
queues,  so  as  to  minimize  the  possibility  of  outdated  or  stale 
measurements.  At  each  node  that  the  BA  packet  traverses 
through,  it  transfers  the  information  that  was  gathered  by  the 
corresponding  FA  packet.  This  information  is  used  to  update 
the  matrices  C  and  1Z  at  the  respective  nodes.  Thus  the  arrival 
of  the  BA  packet  at  the  nodes  triggers  the  iterative  learning 
algorithms.  Of  the  various  learning  algorithms  that  have  been 
proposed  in  the  literature,  we  consider  the  one  proposed  by 
Bean  and  Costa  [2],  In  the  following  sections  of  this  paper, 
we  analyse  this  scheme  for  the  simple  problem  involving 
routing  of  incoming  traffic  between  an  origin-destination  pair 
through  N  parallel  paths. 

Bean  and  Costa  suggest  the  following  scheme  for  the 
learning  algorithms.  Suppose  that  an  FA  packet  measures 
the  delay  A*-0  associated  with  a  walk  ( ).  When 
the  corresponding  BA  packet  arrives  at  node  i  the  delay 
information  is  used  to  update  the  estimate  of  the  mean  delay 
X'j 13  using  the  simple  exponential  estimator 

X)D  :  X)D  +  //!  A'/J  -  A;").  (1) 

where  rj  >  0  is  a  small  constant.  The  mean  delay  estimates 
X1® ,  corresponding  to  the  other  neighbors  to  of  node  i,  are 


left  unchanged. 

Simultaneously,  the  routing  probabilities  at  the  nodes  are 
updated  using  the  scheme: 


j  €  Af(i,  D), 


(2) 


where  /?  is  a  constant  positive  integer.  (3  influences  the 
extent  to  which  outgoing  links  with  lower  delay  estimates  are 
favored  compared  to  the  ones  with  higher  delay  estimates. 

We  can  interpret  the  quantity  as  analogous  to  the 
pheromone  content  on  the  outgoing  link  (i,j).  Equation  (2) 
shows  that  the  outgoing  link  (i,j)  is  more  desirable  when 
XjD,  the  delay  in  routing  through  j ,  is  smaller  (i.e.,  when 
the  pheromone  content  is  higher). 


III.  THE  N  PARALLEL  PATHS  MODEL 

The  model  that  we  consider  pertains  to  the  routing  scenario 
where  arriving  traffic  at  a  single  source  node  S  has  to  be 
routed  to  a  single  destination  node  D.  There  are  N  available 
parallel  paths  between  the  source  and  the  destination  node 
through  which  the  traffic  could  be  routed.  The  network  and 
its  equivalent  queueing  theoretic  model  are  shown  in  Figures 
1  and  2  respectively.  The  queues  represent  the  output  buffers 
(which  we  assume  to  be  infinite)  at  the  source  and  are 
associated  with  the  N  outgoing  links.  We  assume  in  our 
model  that  the  queueing  delays  dominate  the  propagation 
and  the  packet  processing  delays  in  the  N  branches.  These 
additional  delay  components  can  be  incorporated  into  our 
model  with  no  additional  complexity,  but  to  keep  the  dis¬ 
cussion  simple,  we  assume  they  are  negligible.  Two  traffic 
streams,  an  ant  and  a  data  stream,  arrive  at  the  source  node  S. 
At  node  S ,  every  packet  of  the  combined  stream  is  routed 
with  probabilities  (f>i, . . . ,  (f>N  (the  current  values)  towards 
the  queues  Qi,  ■  ■  ■  ,Qn,  respectively.  These  probabilities 
are  updated  dynamically  based  on  running  estimates  of 
the  means  of  the  delays  (waiting  times)  in  the  N  queues. 
Samples  of  the  delays  in  the  N  queues  are  collected  by  the 
ant  packets  (these  are  forward  ant  packets)  as  they  traverse 
through  the  queues.  These  samples  are  then  used  to  construct 
the  running  estimates  of  the  means  of  the  delays  in  the  N 
queues.  We  now  describe  our  model  in  detail. 

We  model  the  arrival  processes  of  ant  and  data  stream 
packets  at  the  source  node  S  as  independent  Poisson  pro¬ 
cesses  of  rates  A  a  and  \p  packets/sec,  respectively.  The 
lengths  of  the  packets  of  the  combined  stream  constitute 
an  i.i.d.  sequence,  which  is  also  statistically  independent  of 
the  packet  arrival  processes.  The  capacity  of  link  i  is  Ct 
bits/sec  (i  =  1, . . . ,  N ).  We  assume  that  the  length  of  an  ant 
packet  is  generally  distributed  with  mean  La  bits,  and  that 
the  length  of  a  data  packet  is  generally  distributed  with  mean 
L I)  bits.  If  we  denote  the  service  times  of  an  ant  and  a  data 
packet  in  queue  Q,  by  the  generic  random  variables  Sf  and 
Sf1,  then  Sf  and  Sf  are  generally  distributed  (according  to 
some  c.d.f.’s,  say  Gf  and  Gf)  with  means  E[Sf]  = 
and  E[Sf]  =  -^P-,  respectively.  The  ant  stream  essentially 


acts  as  a  probing  stream  in  our  system  collecting  samples 
of  delays  while  traversing  through  the  queues  along  with  the 
data  packets.  Thus,  the  packets  of  this  stream  would  be  much 
smaller  in  size  compared  to  the  data  packets. 


Ant  Stream  Capacity  C  , 


Fig.  1.  The  network  with  N  parallel  paths 


Ant  Stream 
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Fig.  2.  N  parallel  paths  :  The  queueing  theoretic  model 


Let  {Aj(m)}  denote  the  sequence  of  delays  experienced 
by  successive  ant  packets  traversing  Qi.  Here  delay  refers 
to  the  total  waiting  time  in  the  system  Qi  (waiting  time  in 
the  queue  plus  packet  service  time).  Also,  let  (5(n)}  denote 
the  sequence  of  successive  arrival  times  of  ant  packets  at  the 
destination  node  D.  Then  the  n-th  arrival  of  an  ant  packet 
at  D  occurs  at  8(n).  Suppose  that  this  ant  packet  has  arrived 
via  Qi.  We  denote  the  decision  variable  for  routing  by  f?(n); 
that  is,  for  *  =  1, . . . ,  N,  we  say  that  the  event  {R(n)  =  i} 
has  occurred  if  the  n-th  ant  packet  that  arrives  at  D  has 
been  routed  via  Qi.  tpi(n)  =  Y^k-i  I{R{k)=i}i  thus,  gives  the 
number  of  ant  packets  that  have  been  routed  via  Qi  among 
a  total  of  n  ant  arrivals  at  destination  D  {I a  is  the  indicator 
random  variable  of  event  A).  Once  the  ant  packet  arrives,  the 
estimate  Xl(n)  of  the  mean  of  the  delay  through  queue  Qi 
is  immediately  updated  using  a  simple  exponential  averaging 
estimator 


Xi(n)  =  Xi(n  -  1)  +  e  (A i(tpi(n))  -  Xx(n  -  1)),  (3) 

0  <  e  <  1,  being  a  small  constant. 

The  delay  estimates  for  the  other  queues  are  left  un¬ 
changed,  i.e.. 


Xj(n)=Xj(n-l),  je{l,...,N},  j±i.  (4) 


In  general,  thus,  the  evolution  of  the  delay  estimates  in  the 
N  queues  can  be  described  by  the  following  set  of  stochastic 
iterative  equations 

Xi(jl)  Xi(jl  1)  T-  £  I{R(n)=i}  (dPi  (tt-)) 

Xi(n  —  1)) ,  i  =  1, . . .  ,N,  (5) 


along  with  a  set  of  initial  conditions  X\  (0)  =  xi , . . . , 

-Ev(O)  =  xN. 

At  time  5(n),  besides  the  delay  estimates,  the  routing 
probabilities  fa(n),i  =  1  , ...,7V,  are  also  updated  simul¬ 
taneously  according  to  the  equations 


fa(n) 


Ef=i  (Xjin))-*’ 


i  =  1, 


,7V, 


(6) 


/?  being  a  constant  positive  integer.  The  initial  values  of  the 
probabilities  are  fa(Q)  =  „ nX ’] — — g -,i  =  1  We 

note  that  fain)  +  •  •  •  +  <j>N(n)  =  1,  for  all  n. 

In  the  above  model  we  consider  only  forward  ants  and  do 
not  incorporate  the  effects  of  backward  ants.  More  precisely, 
we  assume  that  the  estimates  Xj(.)  of  the  means  of  the 
delays  and  the  probabilities  fa(.)  are  updated  as  soon  as 
the  (forward)  ant  packets  arrive  at  the  destination  D,  and 
this  information  is  available  instantaneously  thereafter  at  the 
source  node  S.  We  do  not  consider  thus  the  additional  delay 
involved  as  the  backward  ant  packets  travel  back  carrying  the 
delay  information  to  the  source.  Because  backward  ants  are 
expected  to  travel  back  to  the  source  through  priority  queues, 
the  effects  of  delayed  information  might  not  be  very  signif¬ 
icant,  except  for  large-sized  networks.  On  the  other  hand, 
incorporating  the  effect  of  delays  in  our  model  introduces 
additional  asynchrony,  making  the  problem  harder. 

Another  important  point  to  note  is  that  our  delay  esti¬ 
mation  scheme  (5)  is  a  constant  step  size  scheme.  As  is 
well  known  in  the  literature  on  adaptive  algorithms  (see,  for 
example  [3]),  this  enables  the  scheme  to  adapt  to  (track) 
long  term  changes  in  statistics  of  the  delay  processes.  This  is 
important  for  communication  networks,  because  the  statistics 
of  arrival  processes  at  the  nodes  as  well  as  the  network 
characteristics  typically  change  with  time. 


IV.  ANALYSIS  OF  THE  ALGORITHM 

We  view  the  routing  algorithm,  consisting  of  equations 
(5)  and  (6),  as  a  set  of  discrete  stochastic  iterations  of 
the  type  usually  considered  in  the  literature  on  stochastic 
approximation  methods  [3],  We  provide  below  the  main 
convergence  result  which  states  that,  when  e  is  small  enough, 
the  discrete  iterations  are  closely  tracked  by  a  system  of 
Ordinary  Differential  Equations  (ODEs).  In  the  Appendix, 
we  provide  a  heuristic  analysis  which  enables  us  to  arrive  at 
the  appropriate  ODE. 


A.  The  ODE  Approximation 

An  analysis  of  the  dynamics  of  the  system,  as  given 
by  equations  (5)  and  (6),  is  fairly  complicated.  However, 
when  e  >  0  is  small,  a  time-scale  decomposition  simplifies 


matters  considerably.  The  key  observation  is  that,  when  e 
is  small,  the  delay  estimates  X,  evolve  much  more  slowly 
compared  to  the  waiting  time  (delay)  processes  A ;.  Also, 
because  the  probabilities  fa  are  (memoryless)  functions  of 
the  delay  estimates  Xi:,  they,  too,  evolve  at  the  same  time- 
scale  as  the  delay  estimates.  Consequently,  with  the  vec¬ 
tor  {X\ (n), . . .  ,Xjv(n))  fixed  at  (zi,  ,f.:,  Zn)  (equivalently, 
fa(n),  i  =  1, . . . ,  N,  fixed  at  fa  = 

2^j=l  \z3  ) 

i  —  1,  . . . ,  TV),  the  delay  processes  Aj(.),  i  =  1,  . . . ,  N, 
can  be  considered  as  converged  to  a  stationary  distribution, 
which  depends  on  Also,  when  e  is  small,  the 

evolution  of  the  delay  estimates  can  be  tracked  by  a  system 
of  ODEs.  A  heuristic  analysis  of  the  algorithm,  as  provided 
in  the  Appendix,  shows  that  the  ODE  system  for  our  case  is 
given  by 


dzi(t) 

dt 


Oi  (*))  0  (Di  {z\  (f  ),•••,  zN  (t))  -  z1  (f )) 

E  (zk(t))~0 

3= 1 


dzN{t) 

dt 


( zN(t ))  '3[DN(zi(t)1 . . .  ,zN(t ))  -  zN(t )) 

^  —R  ’ 

E  0 

3=1 

(7) 


with  the  set  of  initial  conditions  Zi(0)  =  xi, . . .  ,Zn{ 0)  = 
xn-  Di(zi, . . .  ,Zn)o  =  1  are  the  mean  waiting 

times  in  the  queues  under  stationarity  (as  seen  by  arriving 
ant  packets)  with  the  delay  estimates  considered  fixed  at 

Zi,...,ZN. 

Formally,  the  ODE  approximation  result  can  be  stated  as 
follows  (see  Benveniste,  Metivier,  and  Priouret  [3]).  For  any 
fixed  e  >  0  and  for  i  =  1 , ,N,  consider  the  piecewise 
constant  interpolation  of  A,(n)  given  by  the  equations  : 
zf(t)  =  Xi(n)  for  t  £  [ne,(n  +  l)e  ),  n  =  0,1,2,..., 
with  the  initial  value  zf( 0)  =  A,(0).  Then  the  processes 
{zf  (t),  t  >  0},  i  =  1, . . . ,  N,  converge  to  the  solution  of  the 
ODE  system  (7)  in  the  following  sense  :  as  e  j  0,  for  any 
0  <  T  <  oo, 

sup  I z-(t)  -  Zi(t) I  0,  (8) 

0  <t<T 

P 

where  — »  denotes  convergence  in  probability. 

In  order  to  obtain  the  evolution  of  the  ODE,  we  need 
to  compute  the  quantities  Di(zi, . . . ,  zjv),  for  our  queueing 
system.  We  recall  that  Di(z i, . . . ,  z jv),  *  =  1, . . . ,  N,  refer  to 
the  means  of  the  waiting  times  as  seen  by  ant  packet  arrivals 
to  the  queues  when  the  delay  estimates  are  considered  fixed 
at  Zi, . . . ,  Zn-  Then  the  routing  probabilities  to  the  N  queues 
are  fa,  =  =  1  We  now  discuss 

2^>j=  1  \Z3) 

how  to  compute  the  quantities  Dt(zi, . . . ,  Zm)  given  our 
assumptions  on  the  statistics  of  the  arrival  processes  and  on 
the  packet  lengths  of  the  arrival  streams. 

Under  such  conditions,  every  incoming  arrival  at  source  S 
from  either  of  the  Poisson  streams  (the  ant  or  the  data  stream) 


is  routed  (independent  of  other  arrivals)  with  probability 
<f>i  towards  queue  (), .  Thus  the  incoming  arrival  process  in 
queue  Qi  (for  each  i)  is  a  superposition  of  two  independent 
Poisson  processes  with  rates  A A<\>i  and  A ofi-  Consequently, 
every  incoming  packet  into  Qi  is,  with  probability  xf+xD  , 
an  ant  packet,  and  with  probability  A  ■  a  data  packet. 
Also,  under  our  assumptions  on  the  statistics  of  the  packet 
lengths  of  the  arrival  streams  and  on  the  arrival  processes, 
all  of  the  queues  evolve  as  independent  M/G/l  queues. 
The  cumulative  incoming  stream  into  Qi  is  Poisson  with 
rate  (Aa  +  A d)4>i,  and  every  incoming  packet’s  service  time 
is  distributed  according  to  the  c.d.f  Gf  with  probability 
Aa+ad  anc*  according  t0  the  c.d.f.  Gf  with  probability 
Aa+Ad  '  further  assume  that  the  queues  are  within  the 
stability  region  of  operation  given  by  the  inequalities  : 


(Aa  +  A  D)(t>iE[Si}  <  1,  i  =  1,  . N,  where  E[Si], 

the  mean  packet  service  time  in  Qi,  is  given  by  = 

AA.E[.SG+Aog[g.;  ].  note  tp,at  average  waiting  time 
in  the  system  as  experienced  by  successive  ant  arrivals 
to  queue  Qi,  is  the  same  as  the  average  waiting  time  in 
Qi  by  the  PASTA  (Poisson  Arrivals  See  Time  Averages) 
property.  Thus,  using  the  Pollaczek-Khinchin  formula  for 
the  average  waiting  time  and  assuming  that  the  queues  are 
stable,  we  finally  obtain  the  expression  for  Di(zi, . . .  ,Zn) 


(i  =  1,...,N): 


Di(z1}...,zN )  =  E[Si] 


(Aa  +  A  p)^E[Sf} 
2(1-(Aa  +  A  D)<t>iE[Si\) 


,  (9) 


where  E[Si]  and 
AA£[gf]+Ap£[Sf 
Aa+Ad 

( Zj 


E[Sf } 
and  E[Sf] 


E,f=i 


rL 

Gi) 


=w- 


given  by  E[Sl]  = 


are 

Aa  +  Ac 


and 


Once  the  expressions  for  Di(zi, . . . ,  Zn)  are  available,  we 
can  numerically  solve  the  ODE  system  (7),  starting  with 
the  initial  conditions  zi(0), . . . ,  Zjv(O).  We  observe  in  our 
simulations  that  if  we  start  the  system  with  initial  conditions 
such  that  we  are  inside  the  stability  region,  the  system  stays 
within  the  stability  region  thereafter. 


B.  Equilibrium  behavior  of  the  routing  algorithm 


We  now  obtain  the  equilibrium  points  of  the  ODE  system 
(7)  which  would,  in  turn,  enable  us  to  obtain  the  equilibrium 
routing  behavior  of  the  system.  In  particular,  we  can  obtain 
the  equilibrium  routing  probabilities  and  the  mean  delays  in 
the  system  under  steady  state  operation  of  the  network.  For  e 
small,  the  steady  state  values  of  the  estimates  of  the  average 
waiting  times  (delays)  in  the  N  queues  are  approximately 
given  by  the  components  of  the  equilibrium  points,  z* ,  of 
the  ODE  system  (7).  The  equilibrium  points  of  the  ODE, 
z* ,  must  satisfy  the  set  of  equations  given  by 


(*!*) 


-0 


N 


E  (*J) 


*\-/3 


=  0, 


j=i 


( ZN ) 


-0 


N 

E 

3= 1 


E  (*J) 


-0 


DN(zl...,z*N)-z*N 


0.  (10) 


The  steady  state  routing  probabilities,  (f>\, . . . ,  f*N,  are  related 
to  the  average  delay  estimates,  z\, . . . ,  z*N,  through  the  equa- 

(z*)~ P 

tions,  f*  =  i ,  ,  i  =  1  Because  we  have 

assumed  that  our  queues  are  in  the  stable  region  of  operation, 
the  steady  state  estimates  of  average  delays  must  be  finite, 
and  so  z*  must  be  finite  for  every  i  =  1 .....  AA  Then 

(z*)~P 

the  steady  state  routing  probabilities,  (f>*  =  , 

E,j=\  \zj ) 

i  =  1,  . . . ,  N,  are  all  strictly  positive.  Equations  (10)  then 
reduce  to  :  z*  =  D;  (zf  . . . ,  z*N),i  =  1, . . . ,  N.  We  also 
notice,  from  equation  (9),  that  for  each  i,  Di(zl, . . . ,  zf,) 
is  a  function  solely  of  </>*,  and  so,  with  a  slight  abuse  of 
notation,  we  denote  it  by  .D,  (</>*).  Then,  utilizing  the  fact 
that  (f> *  =  * ,  ,,-0 ,  we  find  that  the  equilibrium  routing 

z^j=i  zj ' 

probabilities,  <f>f, (f>*N,  must  satisfy  the  following  fixed- 
point  system  of  equations 


(Di^py0 

1=1 


(DN{f*N))-0 

E  (Dj(^)rp 

i= i 


(11) 


Notice  that  </>)'  +  •••  +  <fi*N  =  1. 

An  important  point  to  note  is  that  the  steady  state  proba¬ 
bilities,  4>  1)  •  •  •  )  ^N'  must  not  only  satisfy  the  above  system 
of  equations,  but  must  all  be  strictly  positive  and  satisfy 
the  following  inequalities  which  are  the  stability  conditions 
for  the  system:  (Aa  +  A D)f*E[Si]  <  1,  i  =  1, . . . ,  N.  We 
now  show  that  the  system  of  equations  (11)  are  actually 
the  necessary  and  sufficient  optimality  conditions  for  an 
optimization  problem  involving  the  minimization  of  a  convex 
objective  function  of  (</>i , . . . ,  0jv)  subject  to  the  above 
mentioned  constraints.  We  show  as  a  consequence  that,  if 
there  exists  a  solution  to  the  set  of  equations  (11)  that 
also  satisfies  the  above  mentioned  constraints,  then  such  a 
solution  is  unique. 

Consider  the  optimization  problem 

Minimize  F(</> i, . . . ,  fN)  =  Yn=\  ft  x[Di(x)f  dx, 
subject  to  fi  +  ■  ■  ■  +  (f> jv  =  1, 

0  <  <  a±. 


0  <  (f>N  <  apr, 

Where  °i  =  (A.+A^.!'^1’---^' 

Let  us  denote  by  C  the  set  defined  by  the  constraints  of 

the  above  optimization  problem  -  the  feasible  set.  It  is  easy 
to  see  that  C  is  a  convex  subset  of  RN .  It  is  possible  that 
the  set  C  is  empty  (for  a  given  set  of  values  of  Aa,  A d,  and 
E[Si],i  =  1, . . . , N),  which  means  that  there  are  no  feasible 


solutions  to  the  above  optimization  problem  in  such  a  case. 
We  assume,  in  what  follows,  that  there  exists  at  least  one 
feasible  solution  to  the  above  optimization  problem,  i.e.,  C 
is  non-empty. 

Before  we  attempt  to  solve  the  optimization  problem,  we 
make  certain  natural  assumptions  on  the  delay  functions 
Di(x),i  =  1  We  assume  that  the  functions  D^x) 

are  positive,  real-valued,  differentiable  and  monotonically 
increasing  on  their  domains  of  definition.  This  holds  true  in 
most  cases  of  interest,  because  when  the  routing  probability 
for  an  outgoing  link  increases,  the  amount  of  traffic  flow 
into  that  link  also  increases,  resulting  in  an  increase  of  the 
delay.  The  following  proposition  characterizes  the  optimal 
solutions  </>*  of  the  above  optimization  problem. 

Proposition  1:  Given  the  above  assumptions  on  the  delay 
functions  Di(x),i  =  1 , ...  ,7V,  a  probability  vector  <j>*  is  a 
local  minimum  of  F  over  C  if  and  only  if  cjf  satisfies  the 
set  of  fixed-point  equations  (11).  <j> *  is  then  also  the  unique 
global  minimum  of  F  over  C . 

Proof:  The  Hessian  of  F  is  a  diagonal  matrix  given  by 

V2^!,...,^)  =  diag([A(<Mf“1{A(&)  + 

PtMih)}),  (12) 


where  D'(-)  denotes  the  derivative  of  A(.).  Under  the  above 
assumptions  on  the  Di(x)’ s,  V2F(</i,  •  •  • ,  </tv)  is  positive 
definite  over  C,  and  so  F  is  a  strictly  convex  function  on 
C.  Consequently,  any  local  minimum  of  F  is  also  a  global 
minimum  of  F  over  C\  furthermore,  there  is  atmost  one  such 
global  minimum  [4], 

If  </*  =  (<f>*,  ■  ■  ■  >  4>*n)  Is  a  l°cal  minimum  of  F  over  C, 
we  must  have  (Proposition  2.1.2  of  Bertsekas  [4]), 


^  Q  Tp 


(13) 


Let  us  fix  a  pair  of  indices  i,j,i  ^  j.  Then  choose 
<; t>i  =  (j>*  +  S  and  <j>j  =  </>*-  S,  and  let  ^  i,j. 

Now,  choosing  S  >  0  small  enough  that  the  vector  (j>  = 
•  •  • ,  </jv)  is  also  in  C ,  the  above  condition  becomes 


or,  ^[Dmf^rADMW 


By  a  similar  argument,  we  can  show  that  (j>j[Dj((p*)]^  > 
(f>*  [Di {<f>*)\ ,3 .  Thus,  the  necessary  conditions  for  <j>*  to  be  a 
local  minimum  are 


=  ■■■  =  (j)*N[DN{(l>*N)f  ■ 


conditions.  Then  for  every  other  vector  cf>  £  C,  we  have 
-</*)  =  0-  So,  the  quantity 

Wi  m &  -  $)  =  )  =  o. 

Then,  because  F  is  convex  over  C ,  by  Proposition  2.1.2 
of  Bertsekas  [4],  <f>*  is  a  local  minimum.  □ 

For  our  model,  it  is  easy  to  check  that  the  functions  Di(x), 
as  given  by  (9),  are  positive,  real-valued,  differentiable  and 
monotonically  increasing  on  their  domains  of  definition. 
Thus,  there  is  a  unique  equilibrium  probability  vector  <j>* 
which  satisfies  the  fixed-point  equations  (11). 

V.  SIMULATION  RESULTS  AND  DISCUSSION 

We  describe  in  this  section  an  illustrative  example.  The 
queueing  system  as  described  in  Section  III  has  been  im¬ 
plemented  using  a  discrete  event  simulator.  We  present  here 
results  for  the  case  when  the  number  of  parallel  paths  is 
N  =  3.  The  step  size  e  and  the  parameter  [3  were  set  at  the 
values  0.002  and  1,  respectively. 

The  ant  and  data  traffic  arrival  processes  are  Poisson  with 
rates  =  1  and  Ad  =  1,  respectively.  For  the  ant  packets, 
the  service  times  in  the  three  queues  are  exponential  with 
means  E[Sf)  =  1/3.0,  E[Sf]  =  1/4.0  and  E[S£]  =  1/5.0, 
respectively.  For  the  data  packets  also,  the  service  times  in 
the  three  queues  are  exponential  with  means  ElS®]  =  1/3.0, 
E[S 2  ]  =  1/4.0  and  E[S®}  =  1/5.0,  respectively.  The  initial 
values  of  the  delay  estimates  in  the  three  queues  were  set  at 
Xi(0)  =  0.8,  X2(0)  =  2.8,  and  X3(0)  =  5.6.  Then,  the 
initial  routing  probabilities  are  </>i (0)  =  0.7,  <(>2(0)  =  0-2, 
and  1/3(0)  =  0.1,  which  ensures  that,  initially,  we  are  inside 
the  stability  region  of  the  queueing  sytem.  We  observed  in 
our  simulations  that  as  the  queueing  system  evolved  over 
time,  never  did  the  system  become  unstable. 

Figures  3,  4  and  5  provide  plots  of  the  interpolated 
delay  estimates  (zf(t),i  =  1,2,3)  in  the  three  queues, 
averaged  over  ten  sample  paths,  versus  the  ODE  approxi¬ 
mation,  z\ (t),  22(f)j  Z3(t),  obtained  by  numerically  solving 
(7).  We  see  that  the  theoretical  ODE  tracks  the  simulated 
delay  estimates  fairly  well.  Figures  6,  7,  and  8  provide 
plots  of  routing  probabilities  (f>\ (n),</2(n)>  and  </>3(n)-  The 
routing  probabilities  converge  to  the  equilibrium  values  <j)\  = 
3/12,  (f>2  =  4/12,$}  =  5/12,  which  is  actually  the  unique 
solution  to  equations  (11).  These  probabilities  are  in  the 
reverse  order  as  packet  service  times  in  the  three  queues, 
with  link  3  having  the  highest  equilibrium  probability,  and 
link  1  the  lowest. 


VI.  CONCLUSIONS 


Combining  this  with  the  normalization  condition,  $£  +  •••  + 
(p*N  =  1,  gives  us  the  system  of  equations  (11). 

The  necessary  conditions  above  can  also  be  written  in  the 
form 

dF  dF 

s*r<«--s57W>- 

We  check  that  these  conditions  are  also  sufficient  for  <j>*  to 
be  a  local  minimum.  Suppose  cj>*  £  C  satisfies  the  above 


We  have  provided  convergence  results  for  an  Ant  Routing 
Algorithm  for  a  simple  network  consisting  of  N  parallel 
paths  between  a  source-destination  pair.  We  have  explicitly 
modeled  the  link  delays  using  a  stochastic  queueing  model, 
and  we  have  studied  a  routing  scheme  where  the  routing 
probabilities  are  updated  based  on  estimates  of  path  delays. 
We  have  also  shown  that  the  equilibrium  routing  probabilities 
are  solutions  of  a  fixed-point  system  of  equations,  which  in 


turn,  form  the  necessary  and  sufficient  optimality  conditions 
for  a  convex  optimization  problem.  We  aim  to  extend  the 
analysis  to  the  network  case,  where  multiple  traffic  streams 
with  different  destinations  share  a  network  of  links. 


following  approximate  equations 

Xi  (n  +  M)  w  Xi  (n)  +  Me  (A  (n,  M)  - 
Qi(n,M)X1(nfj, 
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Appendix 

The  heuristic  analysis  provided  below  is  largely  inspired 
by  Benveniste,  Metivier,  and  Priouret  (Section  2.2,  Chapter 
2)  [3],  Let  us  consider  the  mean  delay  estimation  scheme 
given  by  equation  (5).  We  can  write,  for  a  positive  integer 
M, 

M 

Xi(n  +  M)  =  Xi(n)  +  e  E  I{R(n+k)= 1} 

k= 1 

(AiOi  (n  +  k))  —  Xi(n  +  k  —  1)), 


M 

XN(n  +  M)  =  V/v(n)  +  e  A  I\R{n+k)=N } 

fc= l 

(Aj v(Vhv(u  +  k))  -  XN(n  +  k-  1)) . 

(14) 


XN(n  +  M)  «  XN(n)  +  Me(pN{n,M)~ 

QN{n,M)XN{nj).  (15) 


Now  we  tty  to  find  approximations  to  the  quantities 

Pi(n,M )  =  ££iiJ{fi(n+fc)=i}A i(Mn+k))  and  = 

v—'  M  t 

— k=1  for  i  =  1  when  the  values  of 

the  mean  delay  estimates  Xi(. Xn(.)  are  considered 
fixed  at  Xi(n), . . . ,  X]y(n),  and  M  is  large.  Then  the 
routing  probability  vector  (cf>i  (.),...,  <^/v(.)),  too,  can  be 
regarded  as  essentially  constant  in  the  interval  {n, . . . ,  n  + 
M},  because  the  probabilities  are  continuous  functions  of 
the  mean  delay  estimates.  The  routing  probabilities  are 
then  approximately  equal  to  <pi(n)  =  = 

E=i  (Aj(nU 

1,  ...,  N.  Now,  assuming  that  M  is  large  enough  that 
a  law  of  large  numbers  effect  takes  over,  the  average 

y'M  j 

Z"fc=1  {M(n+fc)=’} ,  which  is  the  fraction  of  ant  packets  that 
have  arrived  at  destination  via  Qi  when  the  routing  prob¬ 
abilities  are  cf>i(n ),  can  be  approximated  by  (f>i(n).  With 
the  routing  probabilities  fixed,  the  delay  processes  A,(.), 
can  converge  to  a  stationary  distribution,  the  mean  un¬ 
der  stationarity  being  denoted  by  Di(X\(n), . . . ,  Xpj(n)). 
The  quantities  J(B(”+t)^-'}A‘^‘  "+tl>  can  then  be 

approximated  by  (f>i(n).Di(Xi(n), . . . ,  Xpf(n)).  Note  that 
Di(Xi(n), . . . ,  X]y(n))  are  the  mean  waiting  times  as  seen 
by  ant  packets.  Employing  the  approximations  as  described 
above,  we  notice  from  (15)  that  the  evolution  of  the  vector 
(Xi(n), . . . ,  X]y(n))  resembles  that  of  a  discrete-time  ap¬ 
proximation  to  the  following  ODE  system  when  e  is  small 
enough. 


dz\(t ) 
dt 


C zi(t ))  ^(-Di  {zx  (t zN(t))  -  zi(t,y 

E  {zd^y13 

j= i 


dzN{t) 

dt 


>j^DN(z1(t),...,zN(t))  -  zN(ty 

E  (zk(t))-0 

j= i 

(16) 


with  the  set  of  initial  conditions  Zi(0)  =  xi, . . .  ,Zn{0)  = 

xN. 


If  e  >  0  is  small  enough,  the  vector  (Xi(n), . . . ,  Ajv(n)) 
can  be  assumed  to  have  not  changed  much  in  the  discrete 
interval  {n,n  +  1  +  M },  and  we  can  write  the 
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Fig.  3.  The  ODE  approximation  for  X\  (n) 


Fig.  6.  The  plot  for  </> \  (n) 


Fig.  4.  The  ODE  approximation  for  X2  (n) 


Fig.  7.  The  plot  for  (j> 2(^1) 


Fig.  5.  The  ODE  approximation  for  X3 (n)  Fig.  8.  The  plot  for  </>3(n) 


